#include <iostream>
#include <vector>
#include <list>


#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1D.h>
#include <TCanvas.h>
#include <TVector3.h>

#include "GeomSvc.h"
#include "TriggerRoad.h"
#include "TriggerAnalyzer.h"

using namespace std;

int main(int argc, char* argv[])
{
  GeomSvc* p_geomSvc = GeomSvc::instance();
  p_geomSvc->init("geometry_R962");

  Double_t mass, xF, x1, x2;
  Double_t weight;
  Double_t x0, y0, z0;

  Double_t px1, py1, pz1, p1;
  Double_t px2, py2, pz2, p2;

  TriggerRoad* p_road = new TriggerRoad(); p_road->clear();
  TriggerRoad* m_road = new TriggerRoad(); m_road->clear();

  TFile* dataFile = new TFile(argv[1], "READ");
  TTree* dataTree = (TTree*)dataFile->Get("save");

  dataTree->SetBranchAddress("weight", &weight);
  dataTree->SetBranchAddress("z0", &z0);
  dataTree->SetBranchAddress("mass", &mass);
  dataTree->SetBranchAddress("xF", &xF);

  dataTree->SetBranchAddress("px1", &px1);
  dataTree->SetBranchAddress("py1", &py1);
  dataTree->SetBranchAddress("pz1", &pz1);
  dataTree->SetBranchAddress("px2", &px2);
  dataTree->SetBranchAddress("py2", &py2);
  dataTree->SetBranchAddress("pz2", &pz2);
  
  dataTree->SetBranchAddress("p_road", &p_road);
  dataTree->SetBranchAddress("m_road", &m_road);

  TriggerAnalyzer* triggerAna = new TriggerAnalyzer();
  triggerAna->init("roads_DY.root");

  Double_t cut = atof(argv[3]);
  for(std::list<TriggerRoad>::iterator iter = triggerAna->getRoadList(-1).begin(); iter != triggerAna->getRoadList(-1).end(); ++iter)
    {
      if(iter->ratio() < cut) 
	{
	  iter->disable();
	}
      else
	{
	  iter->enable();
	}
    }

  for(std::list<TriggerRoad>::iterator iter = triggerAna->getRoadList(1).begin(); iter != triggerAna->getRoadList(1).end(); ++iter)
    {
      if(iter->ratio() < cut) 
	{
	  iter->disable();
	}
      else
	{
	  iter->enable();
	}
    }

  Double_t w_target = 0.;
  Double_t w_dump = 0.;
  Double_t w_target_LR_acc = 0;
  Double_t w_dump_LR_acc = 0.;
  Double_t w_target_TB_acc = 0;
  Double_t w_dump_TB_acc = 0.;

  TFile* saveFile = new TFile(argv[2], "recreate");
  TH1D* z_all_target = new TH1D("z_all_target", "z_all_target", 25, -200., 200.);
  TH1D* z_all_dump = new TH1D("z_all_dump", "z_all_dump", 25, -200., 200.);
  TH1D* z_LR_acc_target = new TH1D("z_LR_acc_target", "z_LR_acc_target", 25, -200., 200.);
  TH1D* z_LR_acc_dump = new TH1D("z_LR_acc_dump", "z_LR_acc_dump", 25, -200., 200.);
  TH1D* z_TB_acc_target = new TH1D("z_TB_acc_target", "z_TB_acc_target", 25, -200., 200.);
  TH1D* z_TB_acc_dump = new TH1D("z_TB_acc_dump", "z_TB_acc_dump", 25, -200., 200.);

  TH1D* mass_all_target = new TH1D("mass_all_target", "mass_all_target", 25, 1., 10.);
  TH1D* mass_all_dump = new TH1D("mass_all_dump", "mass_all_dump", 25, 1., 10.);
  TH1D* mass_LR_acc_target = new TH1D("mass_LR_acc_target", "mass_LR_acc_target", 25, 1., 10.);
  TH1D* mass_LR_acc_dump = new TH1D("mass_LR_acc_dump", "mass_LR_acc_dump", 25, 1., 10.);
  TH1D* mass_TB_acc_target = new TH1D("mass_TB_acc_target", "mass_TB_acc_target", 25, 1., 10.);
  TH1D* mass_TB_acc_dump = new TH1D("mass_TB_acc_dump", "mass_TB_acc_dump", 25, 1., 10.);

  TH1D* xF_all_target = new TH1D("xF_all_target", "xF_all_target", 100, -1., 1.);
  TH1D* xF_all_dump = new TH1D("xF_all_dump", "xF_all_dump", 100, -1., 1.);
  TH1D* xF_LR_acc_target = new TH1D("xF_LR_acc_target", "xF_LR_acc_target", 100, -1., 1.);
  TH1D* xF_LR_acc_dump = new TH1D("xF_LR_acc_dump", "xF_LR_acc_dump", 100, -1., 1.);
  TH1D* xF_TB_acc_target = new TH1D("xF_TB_acc_target", "xF_TB_acc_target", 100, -1., 1.);
  TH1D* xF_TB_acc_dump = new TH1D("xF_TB_acc_dump", "xF_TB_acc_dump", 100, -1., 1.);
  TH1D* angle_all_target = new TH1D("angle_all_target", "angle_all_target", 25, -0.05, 0.6);
  TH1D* angle_all_dump = new TH1D("angle_all_dump", "angle_all_dump", 25, -0.05, 0.6);
  TH1D* angle_LR_acc_target = new TH1D("angle_LR_acc_target", "angle_LR_acc_target", 25, -0.05, 0.6);
  TH1D* angle_LR_acc_dump = new TH1D("angle_LR_acc_dump", "angle_LR_acc_dump", 25, -0.05, 0.6);
  TH1D* angle_TB_acc_target = new TH1D("angle_TB_acc_target", "angle_TB_acc_target", 25, -0.05, 0.6);
  TH1D* angle_TB_acc_dump = new TH1D("angle_TB_acc_dump", "angle_TB_acc_dump", 25, -0.05, 0.6);


  for(int i = 0; i < dataTree->GetEntries(); i++)
    {
      dataTree->GetEntry(i);

      TVector3 p_vec(px1, py1, pz1);
      TVector3 m_vec(px2, py2, pz2);
      Double_t angle = p_vec.Angle(m_vec);

      if(z0 > 0) 
	{
	  z_all_dump->Fill(z0, weight);
	  mass_all_dump->Fill(mass, weight);
	  xF_all_dump->Fill(xF, weight);
	  angle_all_dump->Fill(angle, weight);
	
	  w_dump += weight;
	}
      else
	{
	  z_all_target->Fill(z0, weight);
	  mass_all_target->Fill(mass, weight);
	  xF_all_target->Fill(xF, weight);
	  angle_all_target->Fill(angle, weight);
	 
	  w_target += weight;
	}

      cout << p_road->getLR() << "  " << m_road->getLR() << endl;  
      if(triggerAna->acceptEvent(*p_road, *m_road) && p_road->getLR()*m_road->getLR() < 0)/////select road //
	{
          if(z0 > 0)
	    {
    	      z_LR_acc_dump->Fill(z0, weight);    
    	      mass_LR_acc_dump->Fill(mass, weight);
    	      xF_LR_acc_dump->Fill(xF, weight);
    	      angle_LR_acc_dump->Fill(angle, weight);
	      w_dump_LR_acc += weight;
	    }	      
	  else
	    {
	      z_LR_acc_target->Fill(z0, weight);
    	      mass_LR_acc_target->Fill(mass, weight);
    	      xF_LR_acc_target->Fill(xF, weight);
    	      angle_LR_acc_target->Fill(angle, weight);
	      w_target_LR_acc += weight;
	    }
	}///if LR


      if(triggerAna->acceptEvent(*p_road, *m_road) && p_road->getTB()*m_road->getTB() < 0)/////select road //
	{
          if(z0 > 0)
	    {
    	      z_TB_acc_dump->Fill(z0, weight);    
    	      mass_TB_acc_dump->Fill(mass, weight);
    	      xF_TB_acc_dump->Fill(xF, weight);
    	      angle_TB_acc_dump->Fill(angle, weight);
	      w_dump_TB_acc += weight;
	    }	      
	  else
	    {
	      z_TB_acc_target->Fill(z0, weight);
    	      mass_TB_acc_target->Fill(mass, weight);
    	      xF_TB_acc_target->Fill(xF, weight);
    	      angle_TB_acc_target->Fill(angle, weight);
	      w_target_TB_acc += weight;
	    }
	}///if TB



  
      p_road->clear();
      m_road->clear();
    }

  cout << "Dump_LR:   " << w_dump_LR_acc << "  " << w_dump << "  " << w_dump_LR_acc/w_dump << endl;
  cout << "Target_LR: " << w_target_LR_acc << "  " << w_target << "  " << w_target_LR_acc/w_target << endl;
 cout << "Dump_TB:   " << w_dump_TB_acc << "  " << w_dump << "  " << w_dump_TB_acc/w_dump << endl;
  cout << "Target_TB: " << w_target_TB_acc << "  " << w_target << "  " << w_target_TB_acc/w_target << endl;

  cout << "Target/Dump before cut: " << w_target/w_dump <<  endl;
  cout << " ==== Target/Dump after LR cut: " << w_target_LR_acc/w_dump_LR_acc << " === after TB cut: " << w_target_TB_acc/w_dump_TB_acc << endl;
  cout << "Target/All before cut: " << w_target/(w_dump + w_target) << endl;
  cout<<" ==== Target/All after LR cut: " << w_target_LR_acc/(w_dump_LR_acc + w_target_LR_acc) << endl;
  cout<<" ==== Target/All after TB cut: " << w_target_TB_acc/(w_dump_TB_acc + w_target_TB_acc) << endl;
  
  saveFile->cd();

  z_all_target->Write();
  mass_all_target->Write();
  xF_all_target->Write();
  angle_all_target->Write();

  z_all_dump->Write();
  mass_all_dump->Write();
  xF_all_dump->Write();
  angle_all_dump->Write();


  z_LR_acc_target->Write();
  mass_LR_acc_target->Write();
  xF_LR_acc_target->Write();
  angle_LR_acc_target->Write();

  z_LR_acc_dump->Write();
  mass_LR_acc_dump->Write();
  xF_LR_acc_dump->Write();
  angle_LR_acc_dump->Write();

  z_TB_acc_target->Write();
  mass_TB_acc_target->Write();
  xF_TB_acc_target->Write();
  angle_TB_acc_target->Write();

  z_TB_acc_dump->Write();
  mass_TB_acc_dump->Write();
  xF_TB_acc_dump->Write();
  angle_TB_acc_dump->Write();




  saveFile->Close();

  /*
  TCanvas* c1 = new TCanvas();
  c1->cd(); hist_z->Draw();
  c1->SaveAs(argv[2]);
  */

  return 1;
}
